Rare region effects in the contact process on networks 
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Networks and dynamical processes occurring on them have become a paradigmatic representation 
of complex systems. Studying the role of quenched disorder, both intrinsic to nodes and topological, 
is a key challenge. With this in mind, here we analyze the contact process, i.e. the simplest model for 
propagation phenomena, with node-dependent infection rates (i.e. intrinsic quenched disorder) on 
complex networks. We find Griffiths phases and other rare region effects, leading rather generically 
to anomalously slow (algebraic, logarithmic, etc. ) relaxation, on Erdos-Renyi networks. We predict 
similar effects to exist for other topologies as long as a non-vanishing percolation threshold exists. 
More strikingly, we find that Griffiths phases can also emerge -even with constant epidemic rates- 
as a consequence of mere topological heterogeneity. In particular, we find Griffiths phases in finite 
dimensional networks as, for instance, a family of generalized small-world networks. These results 
have a broad spectrum of implications for propagation phenomena and other dynamical processes 
on networks, and are relevant for the analysis of both models and empirical data. 



I. INTRODUCTION 



Complex networks constitute a useful unifying eoncept 
with many interdiseiplinary realizations ranging from the 
World Wide Web and other technological or infrastruc- 
ture networks to genetic, metabolic, ecological, or social 
networks [l| . Many efforts have been devoted to elucidate 
non-trivial topological traits of network architectures; in 
particular, networks with scale- free and/or small- world 
properties received a vast amount of attention. In recent 
years, the research focus shifted to dynamical processes 
occurring on them 0-13 ■ Particularly interesting are 
spreading or transport processes which represent a vast 
variety of propagation phenomena occurring on networks: 
microbial epidemics, computer viruses, rumor spreading, 
or signal propagation in neural nets are some examples. 

As it is well-known in Statistical Physics, the presence 
of quenched disorder usually affects the universal behav- 
ior of phase transitions. Quenched disorder may also 
generate novel phases unheard-of in pure systems (both 
in equilibrium and non-equilibrium situations), as is the 
case of Griffiths phases (GP). These are extended regions 
appearing within the disordered phase and character- 
ized -among some other prominent features- by generic 
anomalously slow dynamics and logarithmic or activated 



scaling at the transition point [5H8|- These effects stem 
from the fact that different rare regions, which can be 
in the ordered phase even if the system is globally dis- 
ordered, emerge in such systems. These regions have a 
broad distribution of relaxation times and the convolu- 
tion of them gives raise generically to slow dynamics. 

Heterogeneity in the intrinsic properties of nodes, 
i.e. quenched disorder, is a natural feature of real net- 
works; node- dependent rates appear in all the examples 
of spreading above (owing to the specificity of the indi- 
vidual immune response, presence of anti-virus software, 
and so forth). Networks with node-dependent intrinsic 
properties (or "fitnesses") have been previously studied 
in the literature @ , but not from the point of view that 
interests us here. Apart from intrinsic node-disorder, 
networks have a structural or topological disorder, since 
nodes are in general topologically not equivalent. One 
can then wonder whether topological disorder by itself 
may induce Griffith phases or similar rare-regions effects. 

The role of intrinsic and topological disorder on the 
overall properties of dynamical processes taking place on 
networks has not been much studied so far. In which 
ways can the node-to-node variability affect the overall 
probability of epidemics to propagate or to become ex- 
tinct? Does disorder in the topology of the network mod- 
ify the basic phenomenology of epidemics? Can novel 
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phases or new qualitative behaviors appear? 

In this paper we tackle the study of quenched disor- 
der -both intrinsic and topological- in dynamical pro- 
cesses on complex networks. For this, we look for rare- 
region effects in the simplest possible epidemic model, 
i.e. the contact process By using different types 

of disorder and various network topologies we report on 
the existence of Griffiths effects, including various non- 
trivial regimes with generic slow decay of activity. In 
particular, we first study a contact process with node- 
dependent rates on Erdos-Renyi networks [T]| and find 
strong rare-region effects below the percolation thresh- 
old. Then we study a contact process with constant rates 
but in disordered topologies: we also report on the exis- 
tence of strong rare region effects in cases in which the 
underlying topology has a finite topological dimension. 
Our conclusions arc expected to go beyond the specific 
examples under consideration, and to apply to different 
models, dynamics and topologies, obeying some minimal 
requirements. We believe that the non-trivial effects of 
disorder uncovered here can shed light on anomalous ef- 
fects observed in many different dynamical problems on 
networks. 

The paper is organized as follows: In Section II we first 
define the models under investigation, the pure contact 
process and the disordered quenched contact process. We 
present a theoretical analysis of the behavior of the latter 
on networks, based on optimal fluctuations arguments, 
and compare it with the results of numerical simulations. 
In Section III wc introduce the Generalized Small- World 
(GSW) networks and the 3-regular random networks and 
illustrate the results of numerical simulations of the con- 
tact process on such topologies. Concluding remarks are 
in Sec. IV. A preliminary account of this work has been 
published in Ref. 



II. QUENCHED CONTACT PROCESS ON 
ERDOS-RENYI AND SCALE-FREE NETWORKS 

A. Contact Process on networks 

Let us consider the pure contact process (CP) [l^, [3 
defined on a generic topology. Each node can be in 
one of two states, either infected (active/1) or healthy 
(inactive/0). An infected node heals at rate fi and, with 
rate A, it infects a randomly selected neighbor. If the se- 
lected neighbor is already infected nothing happens. In 
the following fi will be fixed to 1, with no loss of gener- 
ality. 

A zero-th order homogeneous mean-field equation for 
the average activity density, 

p{t)^-pit) + Mm-pit)), (1) 

predicts an absorbing phase transition where the in- 
fection and healing rates compensate each other: i.e. 



xi^^ = 1 , and a decay p(t) ^ at criticality. A slightly 
more refined calculation (heterogeneous mean-field ap- 
proach [l3]) takes into account the fact that the average 
activity density pk depends on the number of connections 
k (degree) of the corresponding vertex. This again pre- 
dicts a transition at Ac^' = 1 if the underlying network is 
uncorrclated (i.e. vanishing degree-degree correlations). 
These results are expected to be exact for infinite dimen- 
sional lattices as well as for fully connected networks. 
Instead for finitely connected networks the threshold is 
shifted to Ac > 1 (as shown in simulations below). This 
occurs because when activity is low it appears in localized 
regions, and this decreases the effective rate of infection 
(i.e. the probability to choose an occupied nearest neigh- 
bor is larger than for random mixing). This effect can 
be taken into account by using a pair- approximation (as 
described in Appendix A) . In the case of a regular graph 
(with all vertices having degree k) it yields an improved 
estimate of the critical point 



Observe that Ac converges to Ac =1 when fc — > oo 
(i.e. for infinite connectivity, for which simple mean-field 
holds) and diverges at the percolation threshold k = 1 
below which the network becomes fragmented 0, [l^ [l^ 
and, consequently, activity cannot be sustained and the 
phase transition disappears. If the network under con- 
sideration is not regular but has some nontrivial degree 
distribution P{k) it is reasonable to expect that the ex- 
pression 

, (2) _ jk) 
~ (fc)-l' 

where (fe) is the average network connectivity, gives a 
good approximation for the threshold, provided the net- 
work is not very heterogeneous (i.e. that P{k) is nar- 
row). In the case considered in our numerical studies 
below, (fc) = 3 and, hence, the critical point for the pure 

CP is predicted to be around Ac = 3/2, in rather good 
agreement with numerical results. 

B. Quenched Contact Process on networks 

We now consider the Quenched Contact Process 
(QCP) 01 i- e. a Contact Process with quenched disor- 
dered infection rate: a fraction 1 — g of the nodes (type- 
I) take a value A and the remaining fraction q (type-II 
nodes) take a reduced value rA, with < r < 1: 

P(A,)-(l-<z)<5(A, -A)+<z<5(A, -rA) . (4) 

Obviously, for q = and q = 1 the model is non- 
disordered with Xc{q = 1) = Xc{q = 0)/r, while for 
< g < 1 one expects Ac to interpolate between these 
limits. In the general disordered case the density can be 
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expressed as p = {1 — q)pi + qp2, where the sub- index 
refers to the node type and piOie the corresponding den- 
sities. At the homogeneous mean-field level, 



Pi{t) = -pi + (1 - Pi)[\{l - q)pi + rXqp2 



(5) 



A standard linear stability analysis leads 
1 — q{l — r)^^ . As in the pure case, this 



for i = 1 
to, Xliq) 

zero-th order result, multiplied by the factor {k)/{{k) — 
1) to account for correlations, provides a good estimate 
for the threshold in generic networks with narrow degree 
distribution 



Al2) = 



{k) 



1 



l-g(l-r)' 



(6) 



Observe that type-I sites exhibit a percolation transition 
where their intrinsic connectivity, (1 — q){k) — 1 [3. IT5|. 
i.e. at qperc = 1 ^ (^)~^- For larger values of q the net- 
work cannot sustain activity if r = 0: type-I clusters are 
finite and type-II ones do not propagate activity. Hence, 
for r = 0, Eq.(l6]) is valid only for q < qperc, while for 
r > it holds generically. 
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FIG. 1. (Color online) Phase diagram for (k) = 3 and r — 
as a function of the spreading rate A and of the fraction of 
type-I nodes, 1 — q. See main text for a detailed description 
of the different phases. 



where P{s) is the probability of having a rare region of 
size s and t(s) is the decay time of activity in such a 
region. 



Understanding the QCP behavior on 
Erdos-Renyi networks 



Case r = 



The behavior of the QCP model on Erdos-Renyi (ER) 
random networks [llj can be predicted by using, as of- 
ten done in disordered systems (see [Gl-fsl flTj). optimal 
fluctuation arguments. These allow to derive the phase- 
diagram depending on the value of the spreading rate A 
and of the fraction q of nodes with reduced infection rate, 
r\ (see Fig. [T] for r = and Fig. |4] for r > 0). 

In what follows, the theoretical predictions for different 
phases are presented and checked against the results of 
numerical simulations of the QCP on ER networks with 
(fc) = 3 (implying qperc = 2/3), and sizes up to TV = 10''. 
Simulations are performed in the standard way [isl . [loj : 
a list of type-I and type-II occupied nodes is kept and the 
total rates rj and ru are calculated. At each time step 
with probability rj /{rj + rjj) a site of type j is randomly 
selected and it either heals (with probability 1/(1 -I- A^)) 
or infects a single randomly selected neighbor provided 
it was empty (with probability Xj / {1 + Xj)). Time is in- 
creased by l/{rj + rjj) and the procedure is iterated. All 
sites are active initially and the global density of active 
nodes p(t), averaged over many runs, is monitored. 



The basic idea of the optimal fluctuation analysis is 
that the long-time decay of p{t) is controlled by the con- 
volution of different rare regions of type-I sites with dif- 
ferent relaxation time. The overall decay can be written 
as the following convolution integral 



Pit) 



j ds s P{s) exp[-i/T-(s)] 



We start considering the case r = 0. Based on the 
different possible functional forms of P{s) and the clus- 
ter density-decay function in the various regions of the 
phase-diagram the following regimes can be predicted 
(see Fig.l): 

(i) Griffiths phase: X > Xc{qperc) and q > qperc- For 
Q > Qperc the network of type-I nodes is fragmented and 
consist of finite clusters, whose size distribution is given 
by HI: 



P{s) 



1 



2TTP 



(8) 



(7) 



where p is the average number of links per node, which in 
our case is p = (fc)q=o(l~9) for type-I nodes. Within any 
given connected cluster of type-I nodes, let us define pioc 
as the local average number of links per node, and -from 
it- an effective local value of q, qioc = 1 — Pioc/{k)q=o. 
Obviously, for connected type-I clusters, qioc < qperc, 
i.e. they are locally above the percolation threshold and 
hence, provided that A > Xciqperc), they are active rare 
regions, where activity survives until a coherent random 
fiuctuation extinguishes it. The characteristic decay time 
t(s) grows exponentially (Arrhcnius law) with cluster 
size, i.e. t(s) ~ to exp(^(A)s), where to and A{X) do not 
depend on s. Plugging Eq. ^ into Eq. ([7]) and using a 
saddle-point approximation, one obtains p{t) ^ f-^(p^^) ^ 
with 0{p,X) = —{p — 1 — \n{p))/A{X); i.e. there is a 
generic power-Iavir decay with continuously vary- 
ing exponents, that is, a Griffiths phase, emerging as 
the result of ''strong" rare-region effects. It is noteworthy 



that next-to-leading corrections provide logarithmic cor- 
rections to the power-laws. Numerical evidence of this 
GP regime is shown in Fig. (lb) of Ref. [T2|). 

(ii) Right at the percolation threshold, q — qperc-, 
p — > 1 and the distribution of finite clusters, Eq. ([5|, 
becomes a power-law. When plugged into Eq. ([7]) the 
contribution of finite clusters leads to a logarithmic de- 
cay with exponent 1/2, p{t) ^ [ln(i/io)]^^^^, expected 
to hold for any A > \c{%erc)^ for which rare regions 
are active. Observe that a finite fraction of sites be- 
longs to the (infinite) percolation cluster, where activity 
can survive indefinitely, therefore, the actual behavior is 
p{t) ~ c(A) -I- [ln(t/to)]~^^^, where c(A) is a constant, 
meaning a discontinuous phase transition here. Strong 
evidence of such a logarithmic decay behavior is given by 
the inset of Fig. (la) of Ref. [H. 

(iii) For q < qperc there is a giant component of type-I 
nodes which, above the critical point given by Eq.®, is 
able to sustain activity in the steady state. At critical- 
ity, a standard mean-field like contact-process behavior 
{p{t) ~ t~^) is expected and has been numerically veri- 
fied (Fig. (la) in [i2])- On the other hand, in the active 
region, the existence of a percolating cluster implies a 
stationary density p{q,\). Besides this giant component 
some other finite type-I clusters do exist: the relaxation 
in such clusters gives rise to anomalously slow relaxation 
toward p{q, A), analogous to that of regime (i) [l7| . 

(iv) Below the GP (region (i)), at some value of A, and 
for any value q > qperc = 2/3, starts a different region. 
Here, the finite clusters, which were locally supercritical 
in phase (i) become typically sub-critical. Nevertheless, 
the connectedness of finite clusters (as well as the local 
control parameter) varies from cluster to cluster and, al- 
though most of them are locally sub-critical, there still 
may form rare clusters with an over-average mean degree, 
which (themselves or parts of them) are locally supercrit- 
ical at a given A. The extension and the characteristic 
extinction time of these rare regions is unbounded (al- 
though large extinction times are very improbable) and 
this circumstance is sufficient to induce a slowcr-than- 
exponential decay of the global density. The distribution 
of extinction times which are related to the geometry of 
rare regions is extremely difficult to estimate analytically. 
Nevertheless is it expected to decay much more rapidly 
than the power law characteristic of phase (i) and, corre- 
spondingly, the density is expected to decay faster than 
any power of t (but slower than exponentially). In this 
case we speak of weak rare-region effects. 

Numerical results, as shown in Fig. [5] and Fig. [3] indi- 
cate a stretched-cxponential decay, p{t) ^ ^-const-t .^^ji^j^ 
exponents a varying from values close to 1, for small A, 
to very small values (converging to 0) for A approaching 
the GP. Note, that in the limit of vanishing exponent the 
stretched exponential becomes a power-law. 



We can use Eq. ^ the other way around and esti- 
mate the kernel as a function of the resulting convolu- 
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FIG. 2. (Color online) Stretched exponential decay of p vs t 
in region (iv) for = 10^, q = 0.9, r — Q and various values 
of A (which increase from top to bottom curves in the figure) . 
As a reference, straight lines correspond to p{f) ~ exp(— ct^'''^) 
(top) and p{t) ~ exp(— ct^'''^) (bottom). 
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FIG. 3. (Color online) Stretched exponential decay of p vs 
t in region (iv) for A'' — 10^, A = 4.2 and various values of 
q (which decrease from top to bottom curves in the figure). 
As a reference, straight lines correspond to p{t) ~ exp(— t^'''^) 
(bottom) and pit) ~ exp(— f) (top), respectively. The closer 
q to the percolation threshold, the larger the probability to 
have long surviving clusters, and hence the slower the decay. 



tion. Indeed, writing the integral in Eq. ^ in terms 
of extinction times rather than the size s of rare re- 
gions and applying the saddle point approximation we 
obtain that an overall stretched exponential decay for 
p{t) (as measured numerically) implies the asymptotic 
form P[t) ^ cxp{~const ■ t"-/(^^°-^) for the distribution 
of extinction times. 

For values of g < 2/3, above the percolation thresh- 
old, the density decay towards the absorbing state is 
still expected to be dominated by the rare finite clus- 
ters which are present beside the spanning cluster and 
to be of strctched-exponential type, as indeed observed 
numerically (results not shown). 
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Percolation threshold for <K>-3 
NON-PERCOLATING PHASE 1 PERCOLATI NG PHASE ( <K>=3 ) 




FIG. 4. (Color online) Phase diagram for r > 0. In contrast 
to the r = case, the active phase extends all the way to 
q = 1. Otherwise, the phase diagram is very similar to the 
one for r — 0, including a Griffiths phase and a region of weak 
rare region effects. 



(v) For A < X{q ~ 0) and any q, no cluster can be 
super-critical, but still different efFective values of t com- 
pete, giving raise, again, to stretched exponential behav- 
ior. For very small values of A, i.e. deep into the absorb- 
ing phase the decay for all clusters is so fast, that the 
distribution of r becomes very narrow and the decay is 
very close to purely exponential (and therefore the expo- 
nent of the stretched exponential becomes very close to 
unity) . 
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FIG. 5. (Color onhne) Decay of p(t) for q = 0.9 and r = 0.05 
(A'' — 10^). The top plot highlights the presence of a generic 
power-law decay. Values of A are (from top to bottom curves) 
11,10.34,10.2,10,9.5,9,8,7,6,4.5,2.7,1. The bottom plot reveals 
the presence of a stretched exponential regime for A < 4.5. 
Values of A increase from top to bottom curves. The straight 
line corresponds to the exponential decay p{t) = exp(— f). 



2. Case r > 

When nodes of type-II conserve some reduced spread- 
ing capability (r > 0) the features of the phase-diagram 
remain essentially unchanged (see Fig. |4]), except for an 
important qualitative modification: even in the phase 
where typc-I nodes do not percolate a global percolation 
is guaranteed by type-II nodes. One consequence is that 
an active phase exists even when q > qperc for values 
of A larger than \c{q) approximately given by Eq. ([B]). 
The Griffiths phase is limited by such a line above and 
Ac('?perc) below. For smaller values of the spreading rate 
A one still expects a stretched exponential decay and, 
below Xc{q = 0) a pure exponential behavior. In Fig. [5] 
results of numerical simulations for r — 0.05 and q = 0.9 
confirm these expectations. Observe the existence of 
some curvature for the decaying curves in the double log- 
arithmic plot; this is due to the existence of logarithmic 
corrections to scaling. 



Summing up, optimal fluctuation arguments explain 
all numerical findings for both r = and r > 0. Rare 



regions play an important role in almost all phase space, 
giving rise to generic slow decay of activity. 

The nature of the transition between the active phase 
and the Griffiths phase is expected to be of activated 
scaling type, i.e. logarithmic decay is expected (see 
Sect, nil B Indeed, our simulation results suggest acti- 
vated critical behavior but computing accurately the cor- 
responding exponents remains an open challenge. Let us 
remark that recently, a strong disorder renormalization 
group calculation has been performed for Erdos-Renyi 
networks, with the conclusion that a strong-disorder fixed 
point emerges even in this infinite dimensional topology 
[2l|. Comparing numerical results with the theoretical 
predictions in [2l| is left for future work. 

D. QCP behavior on scale-free networks 

It is possible to predict what happens when the QCP 
dynamics takes place on some other, more complex, 
topology. We hypothesize that strong rare-region effects 
(i.e. CP) occur in the fragmented phase of networks with 
a finite percolation threshold: if over-active sites cannot 
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FIG. 6. (Color online) Time decay of the average activity p for 
the QCP on a scale- free network built using the Uncorrelated 
Configuration Model The degree distribution decays as 
P{k) ~ k-"^-^. Solid lines are for A'^ = 10^, dashed lines are 
for = 10*. "Values of A decrease from top to bottom curves. 
Notice that, while for A'^ = Iff* a slow decay occurs at the 
transition, as the system size is increased { N = 10^) the 
typical behavior of the pure systems is observed. 



form rare isolated regions, then rare-region effects do not 
appear. Example of such networks are the ransom ER 
above or structured scale- free networks [l^] . On the other 
hand, if the percolation threshold decreases with system 
size and vanishes in the thermodynamic limit (as is the 
case in standard scale- free networks [11]) only weak rare- 
region effects are expected. 

Numerical evidence of this last assertion is provided 
in Fig. [6l It corresponds to a simulation of the QCP 
on scale-free networks ("generated using the uncorrelated 
configuration model [23[). Apparently, power-law generic 
decay is observed when the network is small {N = 10"*); 
however, as the size is increased (TV = 10^) the transition 
becomes similar to the usual transition for the pure con- 
tact process, with no track of generic power-law decay, 
i.e. no evidence of a Griffiths phase. 



III. EFFECTS OF TOPOLOGICAL DISORDER 
A. General considerations 

In what follows, we assume that the infection rate A.^ 
for site i is identical for all sites and explore the possibility 
of rare region effects induced solely by the topological 
irregularities. It is important to remark that, still, the 
infection rate through any given link is non-homogeneous 
if the degrees of sites are heterogeneous. 

The "local" critical control parameter in any given re- 
gion depends on the connectedness of sites in that region 
(see Eq. ([2])), and therefore heterogeneous networks are 
susceptible to exhibit rare region effects: clusters with 
an over-average local connectedness would have a lower 
local critical control parameter and hence they could be 



locally active even if the system is globally in the absorb- 
ing phase. 

The effects of topological disorder are less clear than 
those of intrinsic node disorder. See, for instance, the 
contradiction between numerical results and the Harris- 
Luck criterion for CP on two-dimensional Voronoi- 
Dclaunay network (the latter predicting topological dis- 
order to be relevant and the former showing the contrary 

M)- 

To shed some light on these problems, let us first con- 
sider the CP on a network with bimodal degree distribu- 
tion, P{k) = p6{k ~ ki) 4- (1 - p)6{k - ki) with fci > fc2, 
where a priori one could expect rare active regions (with 
an over-density of fci-nodes) to exist. However, numeri- 
cally we find just conventional, non-disordered, behavior 
with no evidence of anomalous effects for such networks. 
What is the reason for this apparent contradiction? 

In o?-dimensional lattices, disorder is known to be ir- 
relevant for sufficiently high d, where each node has a 
large number of neighbors and the law of large numbers 
precludes rare regions from existing: each site effectively 
"sees" a well defined mean field, homogenous across the 
system. The topological dimension is an extension of the 
concept of Euclidean dimension to arbitrary graphs. It 
measures how the total number of nodes in a local neigh- 
borhood grows as a function of the topological distance 
from an arbitrary root: N{1) ~ l^ [l^. For small- world 
networks (like the ER graph), where N{1) grows (at least) 
exponentially with I and, consequently, the topological 
dimension is formally infinite. On the other hand, for 
disconnected graphs with no macroscopic (giant) compo- 
nent, like the ER graph below the percolation threshold, 
D = 0. 

For networks with D = oo (as ER graphs above the 
percolation threshold or the bimodal graphs above), the 
number of nodes in a local neighborhood is large. There- 
fore wc conjecture, in analogy with the case of lattices, 
that CPs cannot exist. In such networks, the "locality" 
(i.e. the very existence of local neighborhoods), which 
is a basic component of the phenomenology of rare re- 
gions is broken. This means that the exponentially grow- 
ing neighborhood reduces the possibility of forming well- 
separated rare regions. In other words, the surface of 
these regions is proportional to their volume if _D = oo 
and the number of external links through this large sur- 
face has to be below-average such that the region is iso- 
lated from the rest. As opposed to this, CPs may exist 
for ER below the percolation threshold, where Z? = 
and the finite components of the network are isolated 
from each other. 

In order to get more insight into the effects of topo- 
logical disorder we have studied the CP on several types 
of networks with finite D by numerical simulation. We 
have studied Generalized Small- World (GSW) networks 
[25l - [33j . which consist of a one-dimensional lattice and 
an additional set of long-range edges of arbitrary, un- 
bounded, length. The probability that a pair of sites 
separated by a distance I is connected by an edge decays 
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with / as 



The case s > 2 



P{1) ~ upv 



(9) 



for large I, where A/" is a normalization factor enforcing 
the mean degree to be finite. 

These networks interpolate between the case s = 0, 
which is similar to ER graphs in the sense that long edges 
exist with a uniform (^-independent) probability (note, 
however, that the underlying one-dimensional lattice en- 
sures that the networks is always connected) and the 
quasi-one-dimensional network with certain fraction of 
next-to-nearest neighbor edges corresponding to s = oo. 
In general, P{1) decreasing with the edge length, /, results 
in an overall tendency toward forming clusters of consec- 
utive sites possessing an over-average number of internal 
links, as occurs in the extremal case s ^ oo. In the lat- 
ter model simple considerations predict the existence of 
a GP. Clearly, for s > 2, links have a strong tendency to 
be local; actually the probability that a site belongs to 
a sub-graph that contains many internal links and has 
no external long edge is finite in the limit iV — !■ oo. This 
probability is exponentially small in the sub-graph size, 
suggesting the possibility of strong rare region effects. On 
the other hand, for < s < 2, the number of sub-graphs 
specified above is only sub-linear in the system size, and 
hence they are more likely to become irrelevant in the 
limit N ^ oo. 

In the following two subsections we study two dif- 
ferent families of networks -non-regular and regular 
respectively- within this class of generalized small- world 
networks. 



B. Non-regular generalized small- world networks 

We have studied non-regular GSW networks which 
were proposed in the context of long-range percolation 
[25| . Consider N nodes, labeled 1,2,...,A^ and let us 
define a distance between nodes i and j, I = min(|i — 
j\, N — \i — j\). All pairs of sites at distance I = 1 are 
connected with probability 1, while pairs with I > 1 are 
connected with a probability 



(10) 



obeying Eq. © for large values of I. 

The topological dimension of these graphs has been 
shown to depend on s. If s > 2, the average length of 
edges is finite, consequently, D = 1, whereas for s < 2, 
the average length diverges, implying that the topolog- 
ical dimension is formally D = oo jl^, HO, H^l- In the 
"marginal" case s = 2, D is finite and depends on the pre- 
factor /3 (actually D grows with f3; see below) ^27,. ,30.., 35}. 



^ for example a sub-graph consisting of consecutive sites and all 
next-to-nearest neighbor links 



For s > 2, long-edges are irrelevant, D = 1, and their 
net effect is to induce an over-density of links for some 
nodes, i.e. to induce a form of local quenched disorder. 

For the one-dimensional QCP a strong disorder renor- 
malization group analysis shows that the critical behavior 
is governed by an infinite-randomness fixed po mt ( IRFP) 
where the dynamics is logarithmically slow |36l.l37|. This 
results have been extended to higher dimensions j2ll [ssj . 
In particular, for s pre ading dynamics (i.e. starting from 
a single active site [3^) the survival probability P{t) and 
the number N{t) of active sites which are averaged over 
the initial site behave as 



p{t) ^ Ht/to) 



Nit) ^ Ht/tor. 



(11) 



Instead, initializing the system from fully active state, 
the density of active sites decays as 



p{t) ^ Ht/toY 



(12) 



At the IRFP in one dimension the critical exponents are 
exactly known S = a = [3 — \/5)/2, 77 = VS — 1 [36} . 

On the grounds of universality we expect the same crit- 
ical dynamics as for the one-dimensional QCP with node- 
dependent transition rates [s^ . In order to check this 
conjecture we have investigated the model with s = 3 
and (3^2 with extremely long simulations (2^^ Monte- 
Carlo steps) on lattices of size N = 10^. We measure 
the density decay starting from a fully active state and 
determine the effective decay exponent (by assuming it 
is a power-law and using local slopes in a log-log plot) as 



acsit) 



\n{t/t') 



(13) 



where p(t) and p{t') arc neighboring data points. Other 
effective exponents are measured analogously. Numerical 
simulations in the sub-critical region confirm the presence 
of a Griffiths phase with algebraic dependence on time in 
the range: A = 2.70 - 2.78 (FigE]). 

Results at the transition point are found to be com- 
patible with Eq. (|T2|) with the Id QCP exponent a ~ 
(3 - %/5)/2 sa 0.382 in the critical point at Ac = 2.783(1). 
Indeed, in Fig. |9] the effective exponents (local slopes) 
Soft (t) against time. As a comparison we also plot the ef- 
fective exponents S^g (t) obtained for the one-dimensional 
QCP on lattices of size L = 10^ with bimodal disorder 
distribution (|1]) (r = 0.3, q ~ 0.2; observe that this com- 
parison makes sense since for the QCP the exact relation 
a = S holds, i.e. the rapidity reversal symmetry is not 
broken [l3|). The critical point of the Id QCP is found 
to be located at Ac = 5.24(1), in agreement with [43|. As 
can be seen, the convergence towards a value compatible 
with the expected asymptotic value is very slow (see red 
dot on the ordinate in the inset of Fig. [9l corresponding 
to the t 00 limit) owing to the considerable time scale 
to in Eqs. pmi2p . As a consequence, in case of an acti- 
vated scaling, if neither Ac nor the critical exponents are 
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FIG. 7. (Color online) Top: Local slopes of the density decay 
in networks with s — 3 and /? = 2 for different values of A in 
the Griffiths phase Bottom: Time-dependence of the survival 
probability in the same networks. Numbers shown correspond 
to lines from bottom to top. 
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FIG. 8. (Color online) Top: Time-dependence of the survival 
probability in seed simulations in the network with s = 3 and 
/3 = 2 (numbers shown correspond to lines top bottom) . Bot- 
tom: The corresponding number of active sites in (numbers 
shown correspond to lines from bottom to top). A straight 
line in these plots correspond to the critical point as with 
activated scaling with 5 ~ 0.382 and 77 ~ 1.236. 
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FIG. 9. (Color online) Density decay results in the network 
with s = 3 and P = 2 illustrating the existence of acti- 
vated (logarithmic) scaling at criticality (numbers shown cor- 
respond to lines from top to bottom). Inset: local slopes of 
the same data (lower curves), and local slopes of the survival 
probability in the critical Id QCP simulations. Slow conver- 
gence towards a value compatible with the IRFP value of the 
Id QCP can be observed. 



known a priori the usual method of inspecting the local 
slopes in order to estimate exponents is mostly useless 
and one has to resort to more complex methods [ssj . 

Finally, results for N{t) and P{t) shown in Fig. [5] are 
also compatible with Eqs. ([TT]) with the IRFP exponent 
values. 



2. The marginal case s = 2 

For s = 2 the topological dimension increases con- 
tinuously with p. In Ref. [35|, D{j3) has been es- 
timated for different values of /3. For the values of 
/3 studied in this paper 0.1,0.2,1 the dimensions are 
1.104(3), 1.212(5), 2.35(2), respectively. We have stud- 
ied the CP in these cases, expecting an absorbing phase 
transition at some /3-dependent critical point Ac(/3). It 
is reasonable to assume that the addition of edges to 
a network lowers Ac, i.e. Ac(/3) is monotonically de- 
creasing with /3. Consequently, the model with some 
fixed /3 (0 < /3 < 00) must be in the active phase, if 
A > Ac(0) = 3.297848(22) (the critical point of the one- 
dimensional CP, see [1j1) and must be inactive, if A < 1 
(the critical point of the complete graph described by the 
mean- field equations, see Eq. ^ ). Therefore one can 
predict: I < Ac(/3) < Ac(0), with a possible CP also in 
this range. Note that the range [1, Ac(/3)] where a CP can 
emerge on the inactive side of the critical point shrinks 
upon enlarging /?, making it difficult to numerically ob- 
serve GPs for large values of /3. On the other hand, the 
smaller the value of /3 (and hence the smaller the mean 
degree of nodes) the weaker its influence of the pure sys- 
tem fixed point and the larger times are needed in the 
simulations in order to see the true asymptotic behavior. 
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We have studied the dynamics of CP on these net- 
works via density decay as well as spreading simulations 
from a single active seed for different j3 values around 
the corresponding phase transition points (system-sizes 
L = 10^ - 10^ and up to <max < 10® Monte Carlo steps). 
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simulations for (3 = 0.2 are in accordance with these ex- 
pectations but without knowing the critical exponents it 
is hard to accurately estimate the critical point location. 
We applied a method of Ref. [S^ , based on the assump- 
tion that the leading correction to scaling comes from the 
time scale to in Eqs. pim2p . which is the same for differ- 
ent quantities. At criticality both ln[P(t)] and ln[A^(t)] 
are thus asymptotically proportional to ln[ln(t/to)]. Plot- 
ting ln[7V(t)] against ln[P(t)] the data at the critical point 
must fit to a straight line with a slope —ri/5 (see Fig. [TT|) . 
On the other hand in the GP (actually in all the absorb- 
ing phase) the slope tends to -1-1 whereas in the active 
phase it tends to — oo. This allows us to obtain a rough 
estimate of the critical point: Ac = 2.85(1). 

Unfortunately, the ratio rj/5 varies rather sensitively 
with A around the suggested critical point. The 
data at A = 2.84, 2.85, 2.86 give the ratio estimates 
1.5(1), 1.9(1), 2.6(1), respectively, making it very difficult 
to obtain reliable estimates (for comparison, remind that 
the corresponding ratio of the Id QCP is 1^/5 w 3.236). 



FIG. 10. (Color online) Density decay in the /3 = 0.1 network 
of linear size L = 10® for A = 1.287, 1.29, 1.291, 1.293, 1.295, 
1.297, 1.3, 1.302, 1.303, 1.304, 1.3048, 1.306, 1.307 (from bot- 
tom to top), illustrating the existence of a Griffiths phase. 
The inset shows the corresponding continuously- varying local 
slopes. 

As can be seen in Fig. [10] for /3 = 0.1, the density 
decays algebraically with A-dependent exponents in an 
extended range of A (results of seed simulations for j3 = 
0.2 also support the existence of a GP in a gree ment with 
the results for density decay (see Fig. 3 of |12|')'). 

The calculated effective exponents aeff(t) for /3 = 0.1 
do not level-off for large times (see insets of Fig. [T0|) 
instead a slow drift proportional to l/\n{t) can be ob- 
served. The functional dependence of the effective expo- 
nent can be well fitted by 

a^s{t)^a-a/\n{t) (14) 

suggesting the presence of logarithmic corrections of the 
form: p{t) ~ t^"ln°(i). The possibility of such loga- 
rithmic corrections to power laws in the GP was already 
pointed out in case of the QCP and arises naturally us- 
ing optimal fluctuation arguments (see above). The rela- 
tion a = 5 seems to hold in accordance with the rapidity 
reversal symmetry of the CP [4l[. The phenomenologi- 
cal theory of the GP (see 0) predicts that the number 
of active sites in surviving samples, i.e. N{t)/ P{t) grows 
as a power of hit, implying rj = —5 Q. By extrapolating 
the effective exponents to t — > oo we have found that this 
relation is satisfactorily fulfilled. 

As discussed above a Griffiths phase is usually accom- 
panied by a logarithmically slow -activated dynamics- at 
criticality (sec in Eqs. (|lim2p V We expect this scenario 
to hold for /3 > 0, with possibly modified critical expo- 
nents compared to the case s > 2. Results of numerical 
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FIG. 11. (Color online) The logarithm of the number of par- 
ticles plotted against the logarithm of survival probability for 
different values of A in networks of size N = 10*^ for s = 2 and 
P — 0.2 (numbers shown correspond to lines from bottom to 
top) A straight line in this plot signals the transition point. 

Having an estimate of Ac at our disposal, we turned 
to the estimation of a, S and rj separately. Assuming 
that the survival probability has the asymptotic time de- 
pendence P{t) ~ const X [ln(t/to)]~''; we obtain that the 
effective exponent 6es{t) = — has the following 

dependence on time: 

As can be seen, the deviation from the true value is con- 
siderable whenever In t is not much greater than In to • A 
similar form can be obtained for a and rj as well. We 
have calculated the effective exponents from numerical 
data and fitted the form in Eq. ([T5| to them in the 
domain l/ln(t) < 0.12 where the other corrections are 
expected to be negligible. These numerical data and the 
fitted curves can be seen in Fig. [T^ The extrapolated 
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critical exponents, which can be read off from the inter- 
section with the y-axis express considerable error, due to 
the uncertainty of the location of the critical point even 
though the large simulation efforts applied. 

The exponent 6 is found to increase with /3, leaving 
the value 5{l3 = 0) w 0.382 of QCP, whereas rj is less 
accurate, decreasing from the value = 0) w 1.236 of 
the QCP. 
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FIG. 12. (Color online) Numerically calculated effective expo- 
nents (5eff (i) and Ooff (t) (solid lines) plotted against 1/ Int and 
the curves given in Eq. (|15[l fitted to them (dotted lines) in 
the vicinity of the critical point for P — 0.1 (left) and P — 0.2 
(right) (parameters shown correspond to lines from bottom 
to top). 

We have also studied the evolution of the growing 
cluster in seed simulations by measuring its diameter: 
R{t) = \/ {J2i it)/J2i '^j(O): where the occupation 

number ni{t) is 1(0) if site i at time t is active(inactive), 
li is the length of shortest path between site i and the 
initial site and (•) denotes disorder average over samples 
where the process is surviving up to time t. According to 
the phenomenological theory of the GP, the spread R{t) 
grows as a power of Int. Actually, the critical point of 
the Id QCP strong disorder renormalization group pre- 
dicts that R{t) ^ {In ty/'f' with the exponent ^ = 1/2 
[361 . Indeed, we have found the behavior 



R(t) - (Int) 



1/^0 



(16) 



at criticality with -0 not far from 1/2 for (3 ~ 0.2. 

The picture above changes as /? is further increased. 
Indeed; our numerical observations show that, as argued 
above, the size of the GP shrinks upon enlarging /3. Ac- 
tually, for /3 > 1 we can no longer observe a GP and 
we obtain conventional power-law dynamics at the sug- 
gested critical point rather than activated dynamics. At 
a numerical level, however, we cannot rule out the possi- 
bility of the presence of a GP of vanishing width which 
persists to exist for any finite /3 and is accompanied by 
activated critical dynamics which can be observed only at 
time scales that are well beyond the realm of numerically 
attainable ones. 

With increasing /3, the mean degree increases and it is 
tempting to compare the numerical results with those for 
a related model which can be regarded as the "annealed" 
counterpart of networks with quenched topological irreg- 
ularities. In that model, referred to as CP with Levy- 
flights [12, H^, activation is possible not only at adjacent 
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FIG. 13. (Color online) Left: Density time-decay in MC sim- 
ulations started from the fully active state in networks with 
s = 2 and P — 1. Right: Effective exponents acff(t) for s = 2, 
/3 = 5 in the vicinity of the transition point (parameters shown 
correspond to lines from bottom to top). 



sites but at any site, chosen at any time step ("annealed 
links" ) with a probability decaying with the distance as 



P(r) 



(17) 



Field theoretical analyzes proved that for dimensions 
above dc = 2a, the critical exponents vary continuously 
with cr, while below such a critical dimension they coin- 
cide with mean-field critical exponents (up to logarithmic 
corrections at d = dc) (4^ . In one dimension, this means 
that non-trivial, cr-dependent exponents are expected if 
(T > 1/2, whereas they are stuck to the mean-field values 
if (7 < 1/2. So, it is reasonable to compare the quenched 
model to the annealed one with the same decay expo- 
nent, i.e. s = a + 1. For a = 1 corresponding to s = 2, 
the numerical estimate of the density-decay exponent is 
a = 0.52(3) [1^. This is quite close to the estimated 
critical exponents of the corresponding quenched model 
with large /3, see Fig. [131 At /? = 5 , for which we have 
the most accurate estimates, a = 0.515(20). 

As can be seen from the numerical results, the behavior 
of the quenched model is similar to that of the annealed 
one for large enough /?. This similarity is, however, de- 
teriorates for small /3, which is easy to understand on an 
intuitive level, since for /3 < 1, there is a diverging num- 
ber of backbone hnks {0[N'^~'^]) [H, [13], over which no 
long-range activation can occur, hence the approximation 
by an effective CP with Levy-flight must be inappropri- 
ate. 



3. The case s < 2 

In the case s < 2 we choose the normalization factor 
J\f in Eq. (jH]) such that the mean degree is (fc) = 3 for 
all values of s. Simulations show that the trend observed 
for s = 2 and large values of /3 is continued. Namely, no 
signs of a GP can be found and the critical dynamics are 
of conventional power-law type. Taking into account the 
strong finite-size corrections as well as the possibility of 
logarithmic corrections at s = 3/2 (which corresponds to 
d ^ dc in the annealed model) the results (not shown) for 
s < 3/2 arc compatible with the mean-field value a = 1 
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of the decay exponent in the contact process with Levy- 
flights. As expected, the estimated critical exponents do 
not seem to depend on the mean degree (k) but on the 
index s. For example at s = 1.75, the estimated value is 
a = 0.75(1) (not shown) in agreement with the expected 
value of the decay ex pon ent of the annealed model at this 
point is a « 0.72(5) 



C. 3- regular random networks 

In the networks studied so far the degree of nodes is 
heterogeneous. In the following we consider networks 
with a topological disorder which is even weaker in the 
sense that the degree of all nodes is the same (3 in our 
case). 

In order to keep the topological dimension finite we 
need the probability of long edges to decay according to 
Eq. ([9]) with s = 2. Such 3-regular random networks can 
be constructed in the following way [s^ . Initially, let us 
have a one-dimensional periodic lattice with N vertices, 
all of them are of degree 2. Vertices of degree 2 will be 
briefly called "free vertices". Let us assume that N is 
even and fc is a fixed positive integer. Now, a pair of 
free vertices between which the number of free vertices is 
k—1 (the number of non-free vertices can take any value) 
is selected randomly with uniform probability from the 
set of all such pairs and this pair is then connected by a 
link. That means, for fc = 1, neighboring free vertices are 
connected, for k = 2 next-to-neighboring ones, etc. This 
step, which raises the degree of two free vertices to 3 is 
then iterated until all vertices become of degree 3. (The 
last 2(fc — 1) free vertices are paired in an arbitrary way.) 
Remarkably, as shown in Appendix B, the probability of 
long-distance edges is given by Eq. ([9]) with s = 2 for all k 
and the pre- factor is /3 = fc/2. The topological dimension 
for fc = 1 has been shown to be D(/c = 1) = 2 while for 
fc = 2 the numerical estimate is D{k = 2) = 2.27(2) (H]. 



1. Results 

As expected, simulation results obtained for k = 1,2 
are in line with those obtained for non-regular random 
networks with s = 2 and different values of (3. For fc = 1 
(corresponding to /3 = 0.5 above) we can observe a GP, 
where the density decays algebraically (up to possible log- 
arithmic corrections) with exponents continuously vary- 
ing with A, see Fig. [Ml and activated scaling (not shown) 
at criticality. 

On the other hand, for fc = 2 (corresponding to /3 = 1 
above), finite size effects are stronger because the di- 
ameters of the networks are smaller. The density de- 
cay results are shown in Fig. [151 As can be seen from 
the local slopes, the existence of GP is questionable, in- 
stead a conventional critical phase transition appears at 
Ac = 2.1583(1) with the decay exponent a ~ 0.52(5), 
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FIG. 14. (Color online) Density decay in the k — 1 network 
of size L = 10*^ from A = 2.53, 2.535, 2.537, 2.541, 2.543, 
2.548, 2.554, 2.558, 2.562, 2.57 (from bottom to top). The 
inset shows the corresponding local slopes. 
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FIG. 15. (Color online) Density decay in the k = 2 network 
for L = lO'^ for A = from A = 2.156, 2.157, 2.1575, 2.1581, 
2.1583, 2.1584, 2.1586, 2.1588, 2.1593, 2.1598 (from bottom 
to top). The inset shows the corresponding local slopes. 



which is again close to the corresponding value of the 
annealed model. 



IV. DISCUSSION 

In summary, aimed at studying the effect of disorder on 
propagation phenomena occurring on networks, we have 
investigated the simple contact process on top of differ- 
ent network architectures. First, we have considered a 
quenched contact process, in which the infection rate is 
node-dependent, and we have analyzed it on Erdos-Renyi 
random graphs. As a simple example we have taken a 
bimodal disorder distribution in which nodes infect their 
neighbors either with high or low (even vanishing) prob- 
ability. Localized rare regions, with an over-average in- 
fection rate can emerge in the network when their proba- 
bility is below the network percolation threshold. In such 
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a case strong rare region effects appear. These include 
a Griffiths phase and a rich phase diagram characterized 
by generically slow decay of activity. The main reason 
behind such anomalous behavior is that rare-regions are 
exponentially rare, but being locally active, they survive 
for exponentially large times. The convolution of these 
two effects lead generically to slow decay. Simple "op- 
timal fluctuation arguments" have allowed us to under- 
stand the rich emerging phase diagram and to charac- 
terize analytically the emerging regimes. In particular, 
we distinguish between "strong" rare-region effects ap- 
pearing when the process is locally super-critical, and 
"weak" rare-region effects, occurring when the process is 
predominantly locally sub-critical. 

Similar effects may appear on other topologies as long 
as the percolation threshold is finite. For instance, for 
standard scale free networks in which the percolation 
threshold is known to vanish in the large system-size 
limit, no GP can exist. 

In the second part of the paper we keep the infection 
rate constant at all nodes, and focus the attention on the 
effect of network topological heterogeneity. We conjec- 
ture that, at least for the dynamical processes we have 
studied, Griffiths phases and other rare-region effects can 
appear if the network topological dimension is finite. 
Otherwise, (i.e. for infinite dimensional architectures) 
the very concept of "local neighborhood" does not make 
sense; the frontier of any cluster covers almost completely 
the whole network. We have carefully analyzed different 
generalized small-world networks with finite topological 
dimension: they consist of a one-dimensional lattice, with 
additional long-distance links, which exist with a proba- 
bility that decays with their length / as /3Z~^. For effec- 
tively short ranged links (i.e. s > 2) long-distance edges 
are mostly irrelevant: the topological dimension remains 
D — 1. Their main effect is to create quenched disor- 
der. It is therefore not surprising that the results of our 
computer simulations show results compatible with the 
one-dimensional contact process with quenched disorder. 

For s < 2, however, no Griffiths phase emerges, and 
the critical behavior is of the conventional type with ex- 
ponents close to those of the one-dimensional CP with 
Levy-flights. In the marginal case, s = 2 and for small 
value of /3 we observe a GP and an activated critical be- 
havior with /3-dependent exponents. By increasing /3 the 
width of GP shrinks and for large enough f3 it seems to 
disappear and the critical behavior is found to be con- 
ventional. 

It would be interesting to study if slow-relaxation and 
other rare-region effects appear for dynamical processes 
other than the CP, such as the voter model [l^ or in 
GSW-s built on higher dimensional regular lattices. 

The general aspects of the results obtained in this work 
might be relevant for recent developments in dynamical 
processes on complex networks such as the simple model 
of "working memory" [4611 , social networks with heteroge- 
neous communities |47| or slow relaxation in glassy sys- 
tems [Hi. 
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Appendix A: Pair approximation for the pure 
Contact Process 

Using the notation of [l^ , let us call: p the probability 
to have an occupied site u = p(l, 1) the probability, to 
find a pair of occupied sites, z = p(0, 0) the probability, 
to find a pair of empty sites w = p(l,0) = p(0, 1) the 
probability to have in a pair 1 occupied and 1 empty 
site. Normalization imposes 



u + z + 2w = 1. 



(Al) 



Using the Bayes rule P(l|l) = u/p; P(0|1) = w/p, 
P(1|0) = w/(l - p), and P(0|0) = z/(l - p) for the 
conditional probabilities. Equivalently, 

w - p(i, 0) = (1 - p)p{m = {i-p){i-p- p(o|o)) 

= {l-p)il-p-z/{l-p)) 

= 1- p-l + u + 2w, (A2) 

from where, again, we obtain w = p — u, leaving only 2 
unknowns, say p and u. 

The death rate of any occupied site is fixed to 1, while 
the "infection rate" of a given site with j occupied neigh- 
bors is proportional to \j /k. More specifically, the tran- 
sition rate for a state with a central empty site and j 
occupied neighbors is given by: 



k\ 



k:i\{k-])\ 
kj\{k-j)\ 



[p(i|o)F [p(o|o)]'=-^" p(o) 



(A3) 



and a similar expression for the death rate. Using this, 
it is straightforward to obtain 



{i-pf-^ j-[kj\{k-j)\ 

^-^l^w{w + zf-' - p 



(A4) 



where the combinatorial factor stands for the number of 
ways in which an empty site can be infected by j occu- 
pied neighbors (with j g [1,^])- The factor j/k stems 
from the infection rate. On the other hand, the nega- 
tive term represents death events (notice that the term 
—p could also be obtained by adding up all the possible 
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contributions from all the possible configurations of its 
neighbors). Using that w + z = 1 — p, we obtain the final 
expression: 



p{t) — Xw — p — X{p — u) — p. 
Similarly, for u: 

A ^ p k\ 



(A5) 



1 



j .... — ^ w^'^ 



pk 1 ^ j\{k-j) 



(A6) 



where the extra factor j reflects the fact that every-time 
a site (having j occupied neighbors) becomes occupied 
(resp. empty) , j pairs are created (resp. annihilated). 

For the first term in the r.h.s. we have not found any 
closed expression accounting for the sum, while, for the 
second one, using that u + w = p, it is possible to obtain 
a simplified form: 



uit) 



u^z''-^-ku. (A7) 



Despite of the cumbersome aspect of Ea. (|A7|) . it is pos- 
sible to perform a linear stability analysis of the steady 
state solution of the set of equations Ea.(|A5P and Eq. (jA7|) 
without finding explicitly their analytical solution. Actu- 
ally, from Ea. (jA5[) . the steady state obeys A(p — u) = p, 
and hence, p — u = p/ \, and 



A- 1 



A 



-P- 



(A8) 



Evaluating Ea. (|A7p at linear order in p, only the term 
I = 1 contributes to the series: 

\-kwz^~^ = ku 
k 

\{p - u){l - 2p + uf-^ ^ ku. (A9) 
Plugging here the result above 

p = k^p (AlO) 



from where 

fc(A-l) = A (All) 
indicating that the critical point is located at 

k 



k-l 



(A12) 



Appendix B: Pre-factors of 3-regular graphs 

Here we calculate the pre-factor /3(fc) appearing in 
the asymptotic expression of the probability P{1) for 3- 
regular random graphs. Consider the constructing pro- 
cedure described in the text and assume that the initial 
one-dimensional lattice is infinitely large. When a new 
edge is created, the number of free vertices is reduced by 
2. So when the fraction of free vertices is reduced by an 
infinitesimal amount from c to c — dc, this corresponds to 
the generation of a fraction dc/2 of the long edges. The 
mean length of effective short edges (^) (i.e. distances 
between neighboring free vertices) is 1/c, so it is plausi- 
ble to assume that the probability distribution of ^ has 
the scaling property iTc{C) = ctt{^c) when c ^ Q. The 
length I of a generated new link is the sum of k short 
edges, therefore 



(0 = k/c 



(Bl) 



and its probability Pc{l) has the same scaling property 
as T^dO- We can write for the probability that in the 
network (after finishing the construction procedure) two 
sites in a large distance I are connected by a link: 

P{1)^ / Pc{l)-^- / cP{lc)dc 



2 2 Jo 



1 



Icq 



xP{x)dx ^ ^{x) = (B2) 

where 1/Z ^ cq ^ 1 and we have used Eq. IBll Thus, we 
obtain 



k 

s = 2 and /? = — . 

2 



(B3) 
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